Membrane simulation models from nm to ^m scale 
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Recent developments in lipid membrane models for simulations are reviewed. To reduce compu- 
tational costs, various coarse-grained molecular models have been proposed. Among them, implicit 
solvent (solvent-free) molecular models are relatively more coarse-grained and efficient for simulat- 
ing large bilayer membranes. On a ^m scale, the molecular details are typically negligible and the 
membrane can be described as a continuous curved surface. The theoretical models for fluid and 
elastic membranes with mesh or meshless discretizations are presented. As examples of applications, 
the dynamics of vesicles in flows, vesicle formation, and membrane fusion are presented. 

PACS numbers: 



I. INTRODUCTION 

An amphiphilic molecule, such as a lipid and detergent, 
consists of hydrophilic ('water- loving') and hydropho- 
bic ('water-fearing') parts. In aqueous solution, these 
molecules self-assemble into various structures depending 
on the relative size of the hydrophilic part; spherical or 
cylindrical-like micelles, bilayer membranes, and inverted 
micelles. In particular, the bilayer membrane of phos- 
pholipids is the basic structure of the plasma membrane 
and intracellular compartments of living cells, where the 
membranes are in a fluid phase and lipid molecules can 
diffuse in two-dimensional space. A vesicle (closed mem- 
brane) is considered to be a simple model of cells and 
also has applications as a drug-delivery system. 

Many interesting behaviors of bilayer membranes have 
been studied such as fluid-to-gel phase transition, lateral 
phase separation, protein insertion, the lysis of lipid vesi- 
cles, and membrane fusion. The thickness of a lipid mem- 
brane is 5nm and cells are ~ lOfiin in diameter. Thus, 
the length scale of these phenomena varies from nm to 
/im. To investigate the bilayer structure and molecular 
interactions, information on a nm scale is necessary. For 
example, the length mismatch of Inm between a protein 
and lipid can change the stability and interactions be- 
tween proteins [l[ . On the other hand, the morphologies 
of lipid vesicles and cells on a ^m scale can be understood 
by continuous surface theories 0, H, 3 • 

Various types of membrane models have been proposed 
for the simulations. They are classified in two groups 
depending on whether the bilayer structure is implic- 
itly or explicitly taken into account. In the first group 
[Figs. [Ijd), (e)], the bilayer membrane is described as a 
smoothly-curved mathematical surface 0, H, Q ■ This as- 
sumption is valid on length scales greater than the mem- 
brane thickness of 5nm. Thus, a unit segment represents 
not a lipid molecule, but a membrane patch consisting of 
hundreds to thousands of lipid molecules. The informa- 
tion about the bilayer properties is only reflected in this 
case in the values of the elastic parameters. A mesh dis- 
cretization with a triangular or square mesh is typically 
employed for the simulation. To allow large deformations 
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FIG. 1: (Color online) Various membrane models, (a) all- 
atom model of C12E2. (b) coarse-grained model with ex- 
plicit solvent constructed based on the all-atom simulation in 
(a) [l^- (c) solvent-free molecular model [8l[. (d) meshless 
membrane model [t^]. (e) dynamically-triangulated mem- 
brane model Typical lengths of unit segments are shown 
in parentheses. Figures (a) and (b) were provided by Shinoda. 



of a fluid membrane, remeshing 0] is necessary. In 
particular, a dynamically-triangulated membrane model 
1^, @ is widely used for the membrane with thermal fluc- 
tuations. An alternative model is a meshless membrane, 
in which particles self-assemble into a membrane by po- 
tential interactions. The meshless models are very suit- 
able for studying membrane dynamics accompanied by 
topological changes. 

In the second group [Figs. [Ha)-(c)], amphiphilic 
molecules are modeled by atomistic or coarse-grained 
(CG) molecular models, and the solvent is taken into 
account either explicitly or implicitly. Although com- 
puter technology is rapidly growing, 50 ns dynamics of 
hundreds of lipid molecules are typical levels of recent 
simulations by the all-atom models. Thus, CG models 
have been developed to reduce the computational costs. 
CG molecular models with the Lenard-Jones potential Q 
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FIG. 2: (Color online) Bond flip of triangulated mesh to pro- 
duce membrane fluidity. 



and with dissipative particle dynamics (DPD) [T3| po- 
tential [llj have been widely used for the explicit-solvent 
simulations. Recently, the potential parameters in the 
CG molecular models are tuned by atomistic simulations 
[13, M, [13, [H, [11, see Figs, [^a) and (b). To fur- 
ther reduce the computational costs, larger CG segments 
(two or three segment particles per amphiphilic molecule) 
are employed and the solvent is implicitly represented by 
an effective attractive potential between the hydrophobic 
segments, see Figs. [Ifc). This type of implicit solvent 
model is typically called the solvent-free model. 

There is a length-scale gap between these molecular 
models and the curved-surface models. For the molecu- 
lar models, the size of the unit segments can be gradually 
varied from lA to Inm. For the curved-surface models, 
a unit segment represents > 5nm, and the thermal un- 
dulations of smaller lengths are neglected. 

Recently, the CG molecular simulations and triangu- 
lated membrane were reviewed in Refs. [13, [H, [3 ^.nd 
in Ref. d, respectively. The phase separation in mul- 
ticomponent membranes was reported by Kumar in this 
JSPSJ special issue. In this paper, we review the re- 
cent developments in the solvent-free molecular models 
and curved surface models [Figs, mc)-(e)] for single- 
component membranes and also describes some new re- 
sults. In Sees, [n] and IIIIl the triangulated and meshless 
membrane models arc described, respectively. The area- 
difference elasticity (ADE) model is newly applied to the 
dynamically-triangulated membrane. The membrane dy- 
namics in flows are presented for both models. In Sec. 
rvl the solvent-free molecular models are described. The 
membrane fusion and the formation of polyhedral vesicles 
are presented as applications. 



II. TRIANGULATED MEMBRANE 



A. fluid membrane 



The curvature energy of a fluid membrane is given 
by UMM 



'-iCi + C2~Cof + KC1C2 dA (1) 



where Ci and C2 are the principal curvatures at each 
point in the membrane. The coefficients k and R, are the 



bending rigidity and saddle-splay modulus, respectively. 
The spontaneous curvature Co vanishes when lipids sym- 
metrically distribute in both monolayers of the bilayer. 
Lipid membranes typically have k = 20fcBr, where k^T 
is the thermal energy [^ . For membranes of fixed topol- 
ogy without edges, the integral over the Gaussian cur- 
vature C1C2 is an invariant and their properties arc in- 
dependent of K. Since the critical micelle concentration 
(CMC) of lipids is very low, the number of lipid molecules 
in a membrane is essentially constant over the typical ex- 
perimental time scales. Also, the osmotic pressure gener- 
ated by ions or macromolecules in solution, which cannot 
penetrate the lipid bilayer, keeps the internal volume es- 
sentially constant. Under the constraints of a constant 
volume V and constant surface area A, stomatocytc, dis- 
cocytc, and prolate shapes give the energy minimum of 
Fcv with Co = for < < 0.59, 0.59 < V"* < 0.65, 
and 0.65 ^V*<1, respectively where the reduced 
volume V* = 3y/(47ri?g) with i?o = {A/iny/^. 

In a dynamically-triangulated surface model [^, [H, 
[13, [2^ of vesicles and cells, the membrane is described 
by N vertices which arc connected by bonds (tethers) to 
form a triangular network. The vertices have excluded 
volumes. The curvature energy is discretized using angles 
of the neighboriiig faces d, m, [H or using dual lattices 
of triangulation [1, [i, [23- To model the fluidity of the 
membrane, bonds can be flipped to the diagonal of two 
adjacent triangles using the Monte Carlo (MC) method, 
see Fig. [21 The membrane viscosity p^nh can be varied 
by the frequency of the bond flip [2^, [2^ . This model is 
widely used to simulate a vesicle with spherical topology. 
The membrane with open edges was also studied |24l . 
[25} . Topological changes can be taken into account by 
a mesh reconnection [3,[33|. However, this reconnection 
is a discrete procedure, so that it cannot smoothly treat 
the dynamics. 

Recently, the dynamically-triangulated model were ap- 
plied to the dynamics of a fluid vesicle in flows by com- 
bining with boundary element method s [31 ^ 32 [or multi- 
particle collision (MFC) dynamics [H, [HTsJlsIl . MFC 
is a particle-based simulation method proposed by Male- 
vanets and Kapral (35| . In a simple shear flow, v = jv^x , 
the vesicle shows three types of motions [23, [2^, [3l|, [23] . 
(i) Tank-treading motion: The vesicle has the constant 
inclination angle 9 with respect to the flow direction, and 
instead, the membrane is rotating, (ii) Tumbling: The 
angle 9 rotates, (iii) Swinging: The angle 9 oscillates 
around zero. At a low shear rate 7, the vesicle tran- 
sits from the tank-treading to tumbling with the increas- 
ing viscosity contrast ryin/?/o or membrane viscosity rymb, 
where ryin and rjo are the viscosities of the internal and 
external fluids. At the higher shear rate 7, a swinging 
phase appears between the tank-treading and tumbling 
phases. These dynamics were also observed in experi- 
ments [H, [13 , and are understood by the theory for an 
ellipsoidal vesicle [H, [H or the perturbation theory for 
a quasi-spherical vesicle [1^, [4^, [4l|. We found shear- 
induced shape transitions; elongating transitions from 
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stomatocyte and discocyte to prolate, as well as a shrink- 
ing transition from a tumbling prolate to a tank-treading 
discocyte [H, [1^ ■ In capillary flow, the fluid vesicle also 
transits from discocyte to prolate with the increasing flow 
velocity 



B. area difference elasticity 

In experiments, several morphologies of lipid vesi- 
cles were observed in the addition to the above shapes: 
pear, peal-neckless, and branched starflsh-like shapes 
[Ij 0, SI, To explain them, the area difference elas- 
ticity (ADE) model is widely used [1, H^l- In a vesicle, 
the area of two monolayers of the bilayer membrane are 
different with AA = h§{Ci + C^jdA, where h is the 
distance between two monolayers. Since the flip-flop of 
lipids (traverse motion between monolayers) is very slow, 
the area difference AAq = (A^out — Ain)ao preferred by 
lipids is typically different from AA, where iVout and A'in 
are the number of lipids in the outer and inner monolay- 
ers, respectively, and ao is the area per lipid. In the ADE 
model, the energy of this mismatch AA — AAq is given 

by 

= i^(AA - AA,f = %^(m - m,)^ (2) 

with the averaged curvature m = (l/2_Ro) §{Ci +C2)dA. 
The curvature normalized by to = 47r of a spherical 
vesicle, is denoted by as m* = m/47r. We simu- 
lated vesicles using the dynamically-triangulated model 
for F — Fcv + -Fade with the harmonic constraint po- 
tentials of A and V. Figure [3] shows snapshots of the 
Brownian dynamics (BD) simulation, where the motions 
of the membrane vertices are given by the underdamped 
Langevin equations. 

Vesicles of the ADE model also exhibit shape transi- 
tions in shear flow. A budded vesicle [see the bottom 
snapshots in Fig. |4] shows a tumbling motion at the low 
shear rate 7 = 0.92 [see red lines in Fig. 3], even for 
Vin/vo = 1 and 77,ub = 0, where prolate and oblate vesi- 
cles at any V* only show tank-treading. A similar tum- 
bling motion was reported for a budded two-component 
vesicle At the higher shear rate 7 = 3.68, the vesi- 
cle is found to transit into a prolate shape, which shows 
tank-treading [see blue lines in Fig. 01 . As the shear rate 
decreases, the vesicle transits back to a budded shape. 
The vesicle dynamics of other shapes, such as starfishes, 
are an interesting problem for further studies. 



C. elastic membrane 

The membranes (shells) of synthetic capsules, virus, 
and red blood cells (RBCs), have the shear elasticity ^. 
We call them elastic membranes. The elastic membrane 
can be modeled as a fixed mesh (network). The relative 




FIG. 3: (Color online) Snapsfiots of a fluid vesicle of the ADE 
model at K = fcade = 20kBT and N = 2000. (a) Three- 
armed starfish-like shape at V* = 0.4 and tUq = 2. (b) Two 
connected buds inside a vesicle at V* = 0.8 and mg — —0.5. 
The front quarter of the snapshot in (b) is removed to show 
the inside. 



importance of the bending and shear elasticities is deter- 
mined by the dimcnsionless Foppl-von Karman number 
7 = ERqI K, where E is the two-dimensional Young mod- 
ulus [4^. The elastic membrane with 12 disclinations has 
a spherical shape at 7 < 150 and an icosahedral shape at 
7 > 200. The bending or shear elasticity is dominant for 
the former or latter condition, respectively. The prop- 
erties of viral capsids are well explained by the elastic 
membrane models 0, H^, . 

The models of an RBC membrane have been inten- 
sively studied P, m, MM mi, m, M. in the 
molecular-scale models M: the spectorin network 
attached on the bilayer membrane is modeled by bonds 
with a worm-like chain (WLC) potential. In the con- 
tinuum models, the Skalak model M| is widely used. 
Recently, the comparison with the optical tweezer exper- 
iment [5JI became a standard method to tune the elastic 
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FIG. 4: (Color online) Time developments of asphericity asp 
and inclination angle 6 of a vesicle in simple shear flow at 
V* = 0.75, K = 20fcBr, fcade = 14^7, ml = 2, r,i,,/r,o = 1, 
JVmb = 0, and = 500. The vesicle exhibits shape transitions 
from a budded shape to prolate at 7 = 3.68 and from prolate 
to the budded shape at 7 = 0.92. The asp hericity is the 
degree of deviation from a spherical shape [l3l| . asp = {(Ai — 
Mf + (A2 - A3)" + (A3 - \^f}/{2Rl), where Ai < A2 < A3 
are the eigenvalues of the gyration tensor of the vesicle and 
i?g = Ai + A2 + A3. The intrinsic time unit is r = tiqR^/k. 
The other simulation parameters are described in Ref. [29| . 
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FIG. 5: (Color online) Diameters of an RBC stretched by op- 
tical tweezers. The symbols represent the experimental data 
from Ref. [s^l- The broken and solid lines represent the sim- 
ulation data at fe2 = and 1, respectively 



parameters Isol. ISll . 153 . Isst . Figure [5] shows the experi- 
mental data [5J] and our simulation data. In the simula- 
tion, the RBC membrane is modeled by triangular net- 
works with 578 vertices, which are connected by a bond 
potential ?7bo„d = (fci/2)(r - rof{l + (fc2/2)(r/ro - l)^} 
with M = (\/3/4)fci = 6 X 10-^iV/m and k = 2 x lO^i^J. 



FIG. 6: (Color online) Snapshots of an elastic vesicle in cap- 
illary flow at K = 20fcBr and fiRo/ksT = 110 [sl)- (a) E)is- 
cocyte shape in slow flow, (b) Parachute shape in fast flow. 
The upper-front quarter of the snapshot in (b) is removed to 
show the inside. The arrows represent the flow velocities. 



The nonlinear term with k2 ~ 1 gives better agree- 
ment with the experimental data, but an arbitrariness 
in the choice of the nonlinear function still remains. The 
shape transitions of RBC with a stomatocytc-discocyte- 
cchinocytc sequence can be reproduced by the elastic 
membrane with the ADE model [ssj . 

The dynamics of the RBCs and elastic capsules have 
been simulated by the elastic membrane combined with 
boundary element methods fsol. IstI. IssI . immersed bound- 

61J, 



ary methods [5l|, iM, 60, 61j, MFC [3||, and DFD 53 1. 
Figure [6] shows the elastic vesicles in capillary flow 34 1 . 
At the small flow velocities, the symmetry axis of the dis- 
cocyte is found not to be oriented perpendicular to the 
cylinder axis. With the increasing flow velocity, the elas- 
tic vesicle transits into a parachute-like shape while the 
fluid vesicle transits into a prolate shape. The transition 
velocity linearly depends on the elasticities fj. and k. The 
results are in good agreement with the experimental data 
[6^ . In simple shear flow, elastic capsules (HlllHjIll] and 
RBCs [5^ IHHl transit from tumbling to tank-treading 
with shape oscillation (swinging) as the shear rate in- 
creases. These dynamics are caused by an energy barrier 
for the tank-treading rotation [66l ]. 



III. MESHLESS MEMBRANE 

The meshless models are particle-based methods, 
which do not need bond connections between neighbor- 
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FIG. 7: (Color online) Parameter ka/k^T and e dependence 
of the bending rigidity k and the line tension F [73|. The 
solid (red) lines represent data for ratt/o" = 1-8 and p* = 6. 
The dashed (blue) lines represent data for ratt/o" = 1-9 and 
p* = 8. The solid and dashed lines in (a) and (c) show linear 
fits. 




FIG. 8: (Color online) Sequential snapshots of membrane clo- 
sure simulated by the meshless membrane model with BD at 
N = 4000 ~ lOOOA^o, r-att/cr = 1.8, p* = 6, k^/k^T = 5, and 
e/k^T = 6 [t^]. The membrane stochastically forms one or 
two vesicles. 

ing particles. All of the proposed meshless methods are 
solvent-fee methods, but can be extended to explicit sol- 
vent versions. 

The first meshless model was proposed by Drouffc et 
al. in 1991 (67j . The particles possess an oricntational 
degree of freedom and interact with each other via three 
potentials: a soft-core repulsion, an anisotropic attrac- 
tion, and a hydrophobic multibody interaction. The par- 
ticles self-assemble into a fluid membrane. Very recently, 
the modified models with pairwise potentials were pro- 




initial vesicle of 500 particles. The membrane is in a fiuid state 
with Tatt/o" = 1.8, p* — 6, ka/ksT — 10, and e/ksT = 3. (a) 
t = 0, {h)t = 15/7, (c) t = 18/7, and (d) t = 22/7. 

posed [H, [1^ . They also have oricntational degrees of 
freedom and can form fluid membranes. 

Recently, we proposed an alternative meshless model 
[70| in which particles possess no internal degrees of free- 
dom — unlike the other models. In our model, the par- 
ticles interact with each other via the potential 

C/ = £(C/rop + C/att)+A:aC/a, (3) 

which consists of a repulsive soft-core potential Urep with 
a diameter a, an attractive potential Uutt, and a curva- 
ture potential All three potentials only depend on 
the positions of the particles. Two types of curvature 
potentials [zO] were proposed on the basis of the mov- 
ing least-squares (MLS) method [zll, IZ^- Here, we only 
describe the simpler potential, which drives the particles 
onto a plane: Ua = J2i Q^pKi'i)- The degree of deviation 
from a plane, the aplanarity Upi, is defined as 

a — 9A1A2A3 

(Al + A2 + A3)(AiA2 + A2A3 -f A3A1) 

where Ai, A2, and A3 are eigenvalues of the weighted gy- 
ration tensor, = Y^ji^j - "g)(/3j - PG)wmis{ri,j) 
with a,/3 = x,y,z and a smoothly-truncated Gaussian 
function w^is{r) (70| . In this model, k and F can essen- 
tially be independently varied, see Fig. [7l The bending 
rigidity k linearly increases with ka and is almost in- 
dependent of e. On the other hand, the line tension T 
linearly increases with e and is almost independent of ka 
for ka/k-QT > 10. The saddle-splay modulus is roughly 
estimated as k ~ —k. 

Although the hydrodynamic interactions arc not 
present in the model itself, these interactions can be 
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taken into account by combining with a particle-based 
hydrodynamic method such as MPC [33| and DPD 
[a, US] • Thus, the hydrodynamic interactions can be eas- 
ily switched on or off. Static equilibrium properties can 
be investigated by BD or MC, which requires much less 
computational time. One of the disadvantages of mesh- 
less methods, compared with the mesh models, is that 
the vesicle volume cannot be constrained, since the exact 
volume is difficult to calculate. However, the volume can 
be naturally fixed if an explicit solvent is introduced. 

The effects of the hydrodynamic interactions for the 
self-assembly and dissolution are investigated by compar- 
ing MPC with BD The hydrodynamic interactions 
are found to speed up the dynamics in both cases. The 
particles self-assemble into discoidal clusters, and clus- 
ters aggregate, and then large clusters form vesicles. To 
quantitatively investigate the closing dynamics to a vesi- 
cle, the simulations of the initially flat membranes were 
performed. For a membrane slightly above the critical 
size 4A^0j all of them close into vesicles via a bowl- like 
shape. At N = No, the line tension energy of a flat disk 
equals the bending energy of a spherical vesicle [tJ . For 
much larger membranes of iV 3> 4iVo, however, the mem- 
brane stochastically forms two vesicles via an ^'-shaped 
conformation, see Fig. [51 The large line tension induces a 
negative surface tension, and the resulting buckled shape 
can grow. The closing vesicle has an oblate shape in the 
BD simulations with large N [see the left-bottom snap- 
shot in Fig. [5]. With the hydrodynamic interactions, the 
pressure of the internal fluid makes this closing shape 
more spherical. 

When the line tension is reduced, a vesicle dissolves 
via an opened pore on the membrane with ^yN(t) = 
y/NjO) - ct/2 for both MPC and BD [zl]. A similar 
pore-opening was observed in the lysis of a lipid vesicle 
by detergents [75| . 

In simple shear flow, a spherical vesicle of the meshless 
membrane elongates into a prolate shape with its de- 
creasing volume. Similar elongation is observed for the 
vesicles consisting of surfactants [t^. At higher shears, 
the vesicle ruptures into two pieces, see Fig. [HI The mem- 
brane rupture is a fundamental process for the formation 
of multi- lamellar vesicles in shear flows ITTI. PtsI. [7911 . 



IV. SOLVENT-FREE MOLECULAR MODEL 
A. model 

The simulations of lipid molecules with an explicit sol- 
vent require the calculation for a large number of water 
molecules in addition to lipid molecules. A small patch of 
flat membranes can be simulated with 30 water molecules 
per lipid by the atomistic models [13] ■ However, more wa- 
ter molecules are needed for the simulations of vesicles, 
since the formation of a vesicle (Fig. [5]) requires large 
solvent space to prevent membrane interactions through 
the periodic boundary conditions of the simulation box. 
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FIG. 10: (Color online) Temperature T dependence of diffu- 
sion constant D of amphiphilic molecules, consisting of one 
hydrophilic segment and one hydrophobic segment in the bi- 
layer membrane for various cutoff densities p* . 



The self-assembly of amphiphilic molecules in dilute solu- 
tions requires many more water molecules. The solvent- 
free models are more efficient tools for the simulations 
which require a larger solvent space. A similar solvent- 
free approach is also frequently used in the simulations 
of polymers and proteins. 

The first solvent-free molecular model, which self- 
assembles into a bilayer membrane without an explicit 
solvent, was proposed by us in 2001 [8l|. Later, the mod- 
ified models were proposed for membranes [s^. [ssl. [s^. [ssl. 
[sgI . [stI . [ssI and for micelles Is^ . An amphiphilic molecule 
is modeled as rigid [111,!!!,]!! or flexible M [sElM [13, 
ii, iil chains with single [H, [H, [H, ll, lH,!!!] or 
double [87[ hydrophobic tails. The molecules interact 
with each other with pairwise [H, [H, [H, [H, [13, [H or 
multibody [U, [sgI. [sqI p otentials. In particular, the po- 
tentials in Refs. 



are constructed from the radial 
distribution function of atomistic simulations. 

A common feature of the models is the requirement of 
an attractive potential between hydrophobic segments. 
This attraction mimics the "hydrophobic" interaction 
(hydrophobic segments dislike to contact with water). 
One of the simple ideas to estimate this interaction is that 
the hydration energy is assumed to be proportional to the 
solvent accessible surface area (ASA) JK),, i)l] . Since the 
calculation of the ASA is a numerically time-consuming 
task, an effective potential, which is a function of the 
local density p of hydrophobic segments, were proposed 
for protein [94I and lipid molecules [8l|. At a low den- 
sity p < p* — I, the potential acts as a pairwise attrac- 
tive potential, but at a high density p > p*, this attrac- 
tion vanishes, in which the segments are assumed to be 
completely surrounded by other hydrophobic segments. 
This multibody potential can produce a very fast lat- 
eral diffusion of molecules and a wide fluid-phase range 
(0.1 < fcer/e < 0.9) iSll]- A fluid membrane is also gen- 
erated by pairwise attractive potentials, but their CMC 
is very high, and the membrane often coexists with iso- 
lated monomers. Although the range of the fluid phase 
becomes larger for wider pairwise potentials [85| . their 
ranges are still relatively small. Figure [TO] shows that 




FIG. 11: (Color online) Sliced snapshots of a triangular- 
prism-shaped vesicle at A'^ = 1000, k-e,T je = 0.2 and Cqcf — 
0.11 ■ The red spheres and yellow cylinders represent 
the hydrophilic and hydrophobic segments of amphiphihc 
molecules, respectively. 
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the lateral diffusion on the membrane becomes faster 
with the decreasing p*, where p* = 00 corresponds to 
a pairwise potential. The attractive pairwise potential 
of other solvent-free models can be extended to density- 
dependent potentials in order to obtain faster dynamics. 

The solvent-free molecular models have been applied 



to a variety of phenomena: self-assembly to vesicles 
871. membrane fusion [ol, [o^, [1^, membrane fission 
9y, [q^], the formation of polyhedral vesicles [IS, pore 
formation [H, [oll , the adhesion o f nan oparticles|94. l99l| , 
the fluid-gel phase transition [H, Il00l | , phase separ ation 
of lipids [HI, protein inclusion i n m embrane [osl . Il00l |. 
and DNA- membrane complexes |lOl| . Figure [TT] shows 
a polyhedral vesicle. When the bending rigidity is large 
compared with the vesicle radius, the lined defects at 
the edges of polyhedrons are induced by the mismatch of 
the membrane curvature with the spontaneous curvature 
Co of the monolayer. The number of edges and vertices 
increases with the increasing Co [IJ . Recently, a similar 
defect was p ropo sed as a kink structure in a symmetric 
ripple phase [102 1. 



In early years, an alternative implicit-solvent molecu- 
lar model, in which t he head se gment is constrained on a 
plane, was employed [l03l . llO-j ]- In this model, the mem- 
brane cannot be deformed or self-assembled. However, 
the diffusive motion of molecules is still present on a fixed 
membrane geometry and it may provide a good reference 
state to investigate the effects of the thermal undulations 
of a mem brane . Re cently, a phantom-solvent model was 
proposed [102| . Il05j |. where an explicit but very simple 
solvent is employed. In this model, solvent particles have 
a repulsive interaction with amphiphilic molecules, but 
do not interact with each other. Thus, the solvent is in 
the ideal-gas equation of state, PV = nk-oT. so that its 
pressure is easily controlled. 



FIG. 12: (Color online) Schematic representation of the fu- 
sion pathways obtained by molecular simulations. The green 
(gray) lines represent the hydrophobic boundaries of two 
monolayers. 



B. membrane fusion and fission 

The membrane fusion and fission are key events in 
various intra- and intercellular processes, such as pro- 
tein trafhcking, fertilization, and viral infection. In an 
endocytosis pathway, a small vesicle pinches off from 
the plasma membrane and fuses with a lysosome. Sev- 
eral fusion mechanisms and fusion intermediate struc- 
tures have been propo sed for the fusion of biological 

and lipid i neni b rane s IIO6 \. Among th em, the stalk 

model [H, [Toi, [l03, llOSl. llOiA [uS iSll is widely ac- 
cepted and is qualitatively supported by experimental 
studies. The first intermediate in this model, a stalk, is 
a hourglass-like structure that connects only the outer 
monolayers of the vesicles [Fig. [T^b)]. In the stalk mod- 
els, two types of pore-opening pathways from the stalk 
intermediate were previously proposed. Radial expan- 
sion of the stalk results in contact between the inner 
monolayers inside the stalk, a trans- mono layer contact 
state [Fig. fl^ c)]. In the original model [l08l |. expan- 
sion of the contact area results in a disk-shaped bilayer 
consisting of both inner monolayers, called a hemifusion 
diaphragm [Fig. 112^ 0)]. and a fusion pore is formed in 
the hemifusion diaphragm [Fig. [T2l(b)^(c)^(c)^('f)]. 
The original model is based only on the calculation of 
the bending e nergy of the monolayers in fusion interme- 
diates. Siegel [l09| | also took into account the interstice 
(void) energy of junctions of monolayers. In the model 
modified by Siegel, the trans-monolayer contact results in 
pore formation, rather than expansion of the contact area 
[Fig. [Ill(b)^(c)^(f)]. Although the free energy barrier 
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fixed 




opening at the lined defects |97| . 

Recently, membrane fusion was simulated by various 
molecu lar mode ls such as the bond-fluctuati on l a ttice 
mod el flTsl . [iTl. explicit-solvent CG models pTvl. flTl 
[TTlll2(l. D PP models [ml, [HI, [l2| , and the atomistic 



model |124| . These simulations also show fusion path- 
ways as shown in Fig. 1121 The hcmifusion diaphragm 
was f ound to a l so be formed from the side-pore of the 
stalk [llTl . Ill9l Il23j |. The membrane fi ssion was also 
simulated in two-component vesicles [45l jl25l Il26l |. Al- 
though the dependence on external forces [95| and surface 
tension (T8l . lll6l . [l22l . ll23j | was investigated, conditions to 
determine the pathways are not fully understood. The 
membrane fusion and fission are induced by proteins in 
living cells. Their molecular mechanisms are not well 
understood and have not yet been simulated. 



FIG. 13: (Color online) Sliced snapshots of vesicles under ex- 
ternal forces, (a) Trans-monolayer contact in the membrane 
fusion [93. (b) Formation of a stalk (cylindrical) structure 
in the stretched vesicle [9^. The violet spheres represent 
nanoparticles. 

was estimated to be too high (- 200fcBT) in Ref. [Toot . 
the barrier is reduced by a more qua ntitative estimation 
of the monolayer geometry [lOTl . Illl| . The molecular tilt 
in the monolay ers w as proposed as the mechanism to fill 
the void space |llOl | . Recently, the free energy landscape 
of the fusion w as also st udied using the self-consistent 
field theory [Tsl. ITTl . [TTat . 

The first molecular simulation of membrane fusion was 
performed by a solvent-free model [sl]. Small vesicles 
spontaneously fuse via the modified stalk pathway [Fig. 
[T21(b)^(c)— ^(f)] at low temperature. However, at high 
temperature, the vesicles fuse via a new pathway, where a 
pore opens on the side of the stalk and then the elongated 
stalk bends to form a fusion pore [Fig. [T2l(b)^(d)^(f)]. 
The adhesion of a nanoparticle was found to promote 
this fusion process (s^l- Recently, the leakage between 
the interior and exterior of a vesicle was observed in 
an experiment |114| |. This leakage supports the side- 
pore pathway. The vesicles pinched by particles form 
the trans-monolayer contact [Fig. fT^Ta)] and then fusion 
directly occurs [Fig. [Ill((c)^(f)] or via the hcmifusion 
diaphragm [Fig. [Il(c)->(e)^(f)] 

Membrane fission occurs via pathways opposite of the 
fusion. In a stretched vesicle, a stalk structure [Fig. 
[TST b)] is formed via the trans-monolayer contact [Fig. 
[I21(f)^(c)^(b)] [9^. On the other hand, the particle 
adhesion induces pore opening on the necked membrane 
and the pore expansion leads to the stalk formation [Fig. 
[T2](f)^(d)^(b)] [oil. In large spontaneous curvature 
of the monolayers induces the vesicle fission via pore- 



V. SUMMARY 

We have presented three types of membrane models: 
the triangulated membrane, mcshless membrane, and 
solvent-free molecular models. The first two models are 
constructed for large scale (/^m) and the last model is for 
a molecular scale (nm). To study the bending deforma- 
tion of membranes, the details on a molecular scale are 
typically not necessary, therefore, the mesh and mcsh- 
less models are suitable and much more efficient than 
the molecular models. On the other hand, to study the 
phenomena including the non-bilayer structure, such as 
membrane fusion and protein insertion, the molecular 
scale cannot be neglected. 

When compared with the phenomena in thermal equi- 
librium, the nonequilibrium phenomena are not well un- 
derstood. In this paper, we presented several dynamic 
behaviors on a ^m-scale in flows. On the other hand, 
the dynamics on a molecular scale were much less ex- 
plored while a few simulation studies |l27l Il28| were re- 
ported. The electric field is another interes ting external 
field, which can open pores on a membrane |129l | and in- 
duce shape deformation of a vesicle [13(| . In living cells, 
the biomembrane is in a complex nonequilibrium environ- 
ment. The membrane models presented in this paper are 
powerful tools to study the membrane behaviors under 
nonequilibrium as well as under equilibrium conditions. 
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